Complex trajectories in chaotic dynamical tunneling 
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We develop the semiclassical method of complex trajectories in application to chaotic dynamical 
tunneling. First, we suggest a systematic numerical technique for obtaining complex tunneling 
trajectories by the gradual deformation of the classical ones. This provides a natural classification of 
the tunneling solutions. Second, we present a heuristic procedure for sorting out the least suppressed 
trajectory. As an illustration, we apply our technique to the process of chaotic tunneling in a 
quantum mechanical model with two degrees of freedom. Our analysis reveals rich dynamics of the 
system. At the classical level, there exists an infinite set of unstable solutions forming a fractal 
structure. This structure is inherited by the complex tunneling paths and plays the central role in 
the semiclassical study. The process we consider exhibits the phenomenon of optimal tunneling: the 
suppression exponent of the tunneling probability has a local minimum at a certain energy which is 
thus (locally) the optimal energy for tunneling. We test the proposed method by comparison of the 
semiclassical results with the results of the exact quantum computations and find a good agreement. 

PACS numbers: 03.65.Sq,03.65.Xp,05.45.-a,05.45.Mt 



I. INTRODUCTION 

An intrinsic feature of quantum physics is the existence 
of processes which are forbidden at the classical level. 
The text-book examples of such processes are tunneling 
and over-barrier reflection in one-dimensional quantum 
mechanics" more involved topics include atom ionization 
processes [l|, chemical reactions @, false vacuum decay 
in scalar field theory Q, etc. Generically, one introduces 
a certain parameter g 2 (the Planck constant H in quan- 
tum mechanics, coupling constant in field theory, etc.) 
measuring the magnitude of quantum fluctuations, and 
finds that the probabilities V of classically forbidden pro- 
cesses behave exponentially as g 2 — ► 0, 



V ^ g^Ac-P/v 2 



(1) 



Here F > is the suppression exponent and the depen- 
dence of the pre-exponential factor on g is indicated ex- 
plicitly. In this paper we adopt the term "tunneling" for 
any process forbidden at the classical level. This includes, 
in particular, the cases of dynamical tunneling 0, 
when the exponential suppression of the process is not 
related to the existence of a potential barrier. 

A powerful tool for the study of tunneling at small g 2 
is provided by the semiclassical methods. Exploiting the 
semiclassical approach, one reduces the problem of com- 
puting the tunneling probability to the problem of finding 
the relevant solution to the classical equations of motion 



in the complex domain [l], 0, H[ , where both the time vari- 
able t and dynamical coordinates are taken to be com- 
plex. The suppression exponent F is then related to the 
classical action S calculated along an appropriate contour 
in the complex time plane. The shape of this contour, as 
well as the boundary conditions imposed on the solution 
at t — » ±oo, are dictated by the quantum numbers of the 
initial and final states. In the simplest cases of under- 
barrier motion (one-dimensional tunneling, false vacuum 
decay) the contour runs along the imaginary time axis 1 , 
and the relevant solution is real along this axis. Such 
a trajectory may be identified with the "most proba- 
ble escape path" in the configuration space [H, 0] , which 
gives some understanding of the classically forbidden dy- 
namics. In other cases, however, the passage of the sys- 
tem through the classically forbidden region of the phase 
space cannot be treated separately from the preceding 
and following real-time evolution, and the analysis of 
complex tunneling solutions in the complex time domain 
is needed. This happens, e.g., in the study of chemical re- 
actions with definite initial quantum state of reactants Q 
or in the investigation of the induced tunneling processes 
in field theory [7|, [8|, |£J, [l0( . The semiclassical techniques 
based on genuinely complex classical solutions received 
the common name of the method of complex trajecto- 
ries. 

It should be pointed out that the application of the 
above method might be highly non-trivial. Major com- 
plications are related to the issues of existence and 
uniqueness of tunneling trajectories, which are basically 
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the complex solutions to a certain boundary value prob- 
lem. First, it may occur that the boundary value prob- 
lem at hand does not have any solutions at all (see, e.g., 
Ref. [H|). Second, there may exist many (sometimes 
an infinite number of) solutions. Some of them may 
well be unphysical and should be rejected. The iden- 
tification of physical solutions relies very much on the 
particular properties of the system under consideration; 
presently there are no trustworthy criteria applicable in 
general (see Refs. EE EE EE EE EH for the attempts to 
find such criteria). Even after the unphysical solutions 
are eliminated, the problem remains to identify the so- 
lution^) which yield the minimal suppression exponent 
and therefore correspond to the dominant contribution 
to the tunneling amplitude. 

The above difficulties become particularly pronounced 
in the case of tunneling in chaotic systems. Appearance 
of chaos is generic for non-linear systems with many de- 
grees of freedom; the topics related to tunne ling in the 
presence of chaos were addressed in Refs. EE Il4 jlH H7l . 
El, El SSI S3. The semiclassical analysis of chaotic 
tunneling is hindered by the existence of an infinite num- 
ber of semiclassical solutions which form a fractal set in 
the complex phase space EE HE HH . The direct analysis 
of this set with the purpose of identifying the physically 
relevant solutions becomes an elaborative task in the sys- 
tems with many degrees of freedom [23| . Presently, the 
semiclassical analysis of the chaotic tunneling processes 
is limited to the special cases when t he p has e sp ace of the 
system can be explicitly visualized [l3j, 12CJ, |22| , or when 
the small sub-class of periodic tunneling orbits is consid- 
ered . Development of generic methods of classifying 
the semiclassical solutions is of great importance [231 ]. 

In this paper we present new method of obtaining and 
classifying the tunneling trajectories, which is applicable, 
in particular, in the case of many-dimensional (D > 2) 
chaotic tunneling. Namely, we consider the processes 
which proceed classically at some values of the initial- 
state quantum numbers, and become exponentially sup- 
pressed at other values. The technique presented be- 
low enables one to obtain complex trajectories describ- 
ing tunneling by starting from the real classical solutions 
and changing gradually the quantum numbers of the ini- 
tial state. Our procedure has two advantages. First, 
it is generic and numerically implementable. Second, it 
provides a natural classification of tunneling trajectories 
based on the analysis of their classical progenitors. The 
latter classification suggests a heuristic method for sort- 
ing out the least suppressed tunneling path. 

As an illustration we apply the above method to the 
problem of scattering in quantum mechanical model with 
two degrees of freedom. The process we study is a par- 
ticular exemple of over-barrier reflection. We calculate 
the suppression exponent for the reflection probability. 
As the test of the method we compare the semiclassical 
results with the "exact" suppression exponent. The lat- 
ter is extracted from the exact wave function obtained 
by solving numerically the Schrodinger equation. The 



results of the two calculations are in good agreement. 

The system under consideration has two distinctive 
features, which are inherent in two wide (intersecting but 
not necessarily identical) classes of tunneling problems. 
We believe that our model is a generic representative of 
both of these classes. 

The first feature is chaoticity. In our model chaos man- 
ifests itself at the classical level as follows. Consider the 
set of initial data giving rise to the classical reflected 
trajectories. We will see that this set falls into an infi- 
nite number of disconnected domains. The boundaries 
of these domains correspond to trajectories which do not 
escape into the final asymptotic region as time goes on, 
but get trapped in the interaction region. The latter tra- 
jectories are unstable: small deviations from them lead 
to either reflected or transmitted solutions. Increasing 
the resolution of initial data reveals that the set of initial 
data corresponding to the trapped trajectories forms a 
fractal similar to the Cantor set. This is the hallmark of 
the so-called irregular (or chaotic) scattering [24|. 

The chaoticity of the classical dynamics has profound 
consequences for the tunneling process. We will show 
that complex trajectories relevant for over-barrier reflec- 
tion are all trapped in the interaction region and thus are 
unstable. Moreover, they inherit the fractal structure of 
the classical trapped solutions which means, in particu- 
lar, that their number is infinite. The complex trajectory 
which contributes most to the tunneling amplitude is a 
descendant of a certain unstable classical solution lying 
on the boundary of the above fractal set. 

The chaoticity of the process manifests itself in the 
exact quantum computations as well. We find that 
the quantum probability of tunneling, instead of being 
smooth function of energy, exhibits large irregular oscil- 
lations. Similar dependence of the tunneling amplitude 
on the parameters of chaotic systems was reported pre- 
viously in Refs. [13 EE HI- At nrst glance, this be- 
havior contradicts to the semiclassical formula ([TJ. One 
observes, however, that the oscillation period scales like 
g 2 when g 2 — > 0, so that the oscillations become indis- 
cernible in the semiclassical limit. In order to extract 
the semiclassical suppression exponent, one smears the 
tunneling probability over several periods. The smeared 
probability does obey the scaling law |T]). 

The second feature of our system is as follows. We 
observe that the process under consideration is classi- 
cally forbidden, and hence exponentially suppressed, at 
arbitrary high energies. Our interest in this property is 
motivated by the studies of similar processes in quan- 
tum field theory. As a matter of fact, the exponential 
suppression at all energies is generic for the field theo- 
retical processes involving creation of some classical ob- 
ject (soliton, bubble of new phase or vacuum configura- 
tion with different topology) in a collision of two highly 
energetic quantum particles EE 
over, it has been shown recently 



25, |2fl |22l|2S|. More- 
101 ] that the method of 



complex trajectories predicts the suppression exponent 
F of the above processes to attain its minimum at a cer- 
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tain "optimal" energy E a , above which F stays constant 
(this behavior of the suppre ssion exponent was conjec- 
tured earlier in Refs. [26|, [27]]). The optimal value F(E ) 
is determined by the complex-valued classical solutions 
with particular properties, see Refs [lOj; these solutions 
are called "real-time instantons" . 

One might question the applicability of the method 
of complex trajectories for the description of the above 
phenomenon of optimal tunneling. Indeed, the proper- 
ties of processes at E » E Q are in many respects different 
from the well-known case of tunneling through a poten- 
tial barrier. In this paper we provide the evidence that 
the semiclassical method is applicable for the description 
of dynamical tunneling independently of how high the en- 
ergy of the process is, or whether there exists a potential 
barrier at all. 

The model of this paper provides a particular exam- 
ple of quantum mechanical system exhibiting the phe- 
nomenon of optimal tunneling. Namely, the suppression 
exponent F(E) depends on energy E in non-monotonic 
way, attaining a local minimum at some energy E a . We 
will find that the minimal value of the suppression ex- 
ponent is indeed given by the method of real-time in- 
stantons. At higher energies the function F(E) grows to 
infinity 2 . 

The paper is organized as follows. In Sec. [TTJ we de- 
scribe the model under consideration and introduce no- 
tations. In Sec. IIIII the classical dynamics of the system 
is analyzed. The semiclassical study of the classically 
forbidden reflections is performed in Sec. IIVI In Sec. [V] 
we present the results of the numerical integration of the 
Schrodinger equation and discuss their comparison with 
the semiclassical results. Sec. I VII contains summary and 
conclusions. 



II. SETUP 

Throughout the paper we illustrate our technique on 
a toy model describing the evolution of a quantum par- 
ticle in two-dimensional harmonic waveguide. Namely, 
we consider the case when the motion of the particle is 
confined to the vicinity of a certain line Y — ^a(gX) by 
the quadratic potential (X and Y stand for the Cartesian 
coordinates of the particle). The equipotential contour 
U(X, Y) = £ is shown in Fig. [TJ the Hamiltonian is 



Y 



H 



p2 



p2 



2m 



u{xx) , 



(2) 



2 This is different from the case of collision— induced tunneling in 
field theory, where the suppression exponent stays constant at 
energies higher than E . The reason is that in the latter case at 
E > E another tunneling mechanism, which is specific for the 
field theoretical setup, comes into play. Namely, above the opti- 
mal point the energy excess (E — E a ) is released by the emission 
of a few hard quantum particles, so that the tunneling transition 
effectively occurs at the optimal energy [ToL |2^|. 
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FIG. 1: The equipotential contour U(X,Y) — £ for the 
waveguide model @ and the directions of the incoming, re- 
flected and transmitted fluxes of particles. 



where 



U(X,Y) 



a(x) = a Q e x "'I 2 



Y 



-a{gX) 



(3) 



In the asymptotic regions X — > ±oo, the variables of 
the model separate, and the motion of the particle 
becomes trivial: oscillations in the Y"-direction are ac- 
companied by the translatory motion along the X-axis. 
The wriggle around the point X introduces nonlinear 
coupling between the degrees of freedom; we refer to this 
part of the configuration space as the interaction region. 
In what follows we use the system of units where 



1 



The rescaling 



X = x/g, P = pg 
brings the Hamiltonian ([2]) into the form 



H 



— (y-a(x)Y 



(4) 



(5) 



(6) 



Equation ([5]) implies that the semiclassical regime in our 
model occurs for g 2 <C 1. Apart from g 2 , the only free 
parameter of the model ([6]) is an, which we set 



a = 0.8 



(7) 



This choice will be explained in Sec. IIIII To avoid con- 
fusion, we remark that the rescaled coordinates x, y will 
be used for the semiclassical analysis of Sees. IIIII IIVI and 
App. [5] while the original ones (X and Y) are exploited 
in the quantum computations of Sec. [Viand App. iBl 

The process we want to investigate is the backward re- 
flection of a particle coming from the right, see Fig. Q] 
Note that though we will refer to this process as over- 
barrier reflection, there is actually no potential barrier in 
our system: the minimum of the potential is zero in any 
transverse section of the waveguide. The incoming quan- 
tum state \£,Af) is completely determined by the total 
energy £ and occupation number M of the Y~-oscillator. 
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The quantity of interest is the total reflection coefficient 
for this state, 



1 



V(S,Af)= lim 

tf—ti^ + OO If — ti 



X 



x£|(/|e-^-^,A0| 
/ 



(8) 



where |/) stands for the basis of reflected waves (P/ > 0) 
supported in the right asymptotic region, and the proper 
normalization of the incoming state has been chosen, 



(£,Af\£',Af') = 2wS(£ - £')<W 



(9) 



At some values of the initial-state parameters £, Af the 
reflection process is classically forbidden. [In particular, 
one expects a particle with high enough translatory mo- 
mentum | Pi | to pass classically to the other side of the 
waveguide ending up in the asymptotic region X — > — oo.] 
In this case the reflection coefficient © is expected to 
obey the semiclassical scaling law ([TJ at g 2 — > 0. Below 
we concentrate on the calculation of the leading suppres- 
sion exponent F as function of £ , Af. 

We do not pretend to describe any concrete experi- 
mental situation by the Hamiltonian @ ; it is chosen as a 
convenient testing ground of our semiclassical technique. 
The advantage of the scattering setup is an unambiguous 
determination of the initial and final states of the tunnel- 
ing process, the latter determination being problematic 
in the case of bounded motion Q . Moreover, the simple 
form of our model makes it tractable both semiclassically 
and by the exact quantum mechanical methods. Note 
that a system similar to © was considered in Ref. [29| . 



III. CLASSICAL REFLECTIONS 

Let us start by considering the classical dynamics of 
the reflection process. As we will see in Sec. HVl the 
analysis of the classical dynamics is crucial for under- 
standing the classically forbidden reflections. The results 
of the present Section will enable us to identify the sets 
of the initial data (£, Af) which correspond to the clas- 
sically allowed/forbidden reflections; for brevity we refer 
to these sets as classically allowed/classically forbidden 
regions. Note that these are the regions in the initial 
data plane (£, Af); they should not be confused with the 
classically allowed/forbidden regions in the configuration 
space. The latter term is not used in this paper. 

One notes that the action functional of the model © 
has the form 



S = S/g 2 , 



(10) 



where 



5' 



x 2 + y 2 1 
2 2 l 



(y - a(x)) 2 dt (11) 



does not contain the parameter g at all. Hence, g drops 
out of the classical equations of motion. While studying 
the classical dynamics it is convenient to forget about g 
and consider the classical system defined by the rescaled 
action (fTTj) . 

In contrast to quantum mechanics, where two quantum 
numbers £ , Af determine the initial state completely, the 
classical evolution is specified by four initial conditions. 
One of these is physically irrelevant. It corresponds to 
the x coordinate at the initial moment t = ti, and can be 
absorbed by the appropriate shift of ti. Note that, still, 
one should be careful to choose x{tf) far from the interac- 
tion region; in numerical calculations of this Section the 
value 



x{U) = 10 



(12a) 



is used. To keep contact with the quantum mechanical 
formulation we choose the other two initial data to be 
the total classical energy E and the "classical occupation 
number" N, which is equal 3 to the initial classical energy 
of the transverse oscillations. This determines the initial 
velocity of the particle along the waveguide, 



x{U) = -y/2{E-N) 



(12b) 



It is worth noting that Eq. (| 10(1 implies the following rela- 
tion between the classical parameters and their quantum 
counterparts, 



E = g 2 £ , N = g 2 Af 



(13) 



The last initial condition is the initial phase 0o of y- 
oscillator. It parametrizes the initial position and veloc- 
ity of the particle in the transverse section of the waveg- 
uide: 



y(ti) = V2N cos 0o , 
y(U) = -%/2iVsin0 o 



(14a) 
(14b) 



Clearly, a given point of the plane (E, N) belongs to the 
classically allowed region if classical reflection is possible 
for some value(s) of 4>q. Otherwise, we say that this point 
lies in the classically forbidden region. 

At first glance it may seem that classical reflections are 
impossible for any values of E, N as there is no poten- 
tial barrier to prevent the classical particle from going 
into the left asymptotic region. Let us make sure that 
this is not the case by considering the evolution at small 
E. In this regime the particle moves slowly along the 
axis y = a(x) of the waveguide performing small and 
(relatively) rapid oscillations in the orthogonal direction. 
The frequency u>± of the latter oscillations is determined 
by the curvature of the transverse section of the poten- 
tial and thus depends on the position x of the particle. 



Recall that in our units i 



■5 



It is straightforward to find that coj_(x) = y/l + (a'(x)) 2 . 
The energy of the orthogonal oscillations E±_ divided by 
their frequency is an adiabatic invariant, 



E\ 



LUj_(x) 



= const = N 



(15) 



On the other hand, the conservation of total energy 
yields, 



E 



+ iV v /l + (a'(x)) 2 



(16) 



where vn is the projection of the particle velocity onto 
the axis of the waveguide. Thus, the adiabatic motion in 
our waveguide is governed by one-dimensional dynamics 
in the effective potential 



U(x) = Nyfl- 



-.2 «,2 £ 



(17) 



where the explicit expression ([3]) for the function a(x) has 
been used. This picture is valid as long as wi/w} <C 1, 
which is satisfied for E,N <C 1. 

The potential (fTT)) has the form of two symmetric 
humps with the maxima U max = Ny/l + a^e -1 situated 
at x = ±1, see Fig [2] Any particle coming from the right 




FIG. 2: The effective potential for the motion in the adiabatic 
regime. The particle is reflected at E < Umax (case a), and 
is transmitted through the waveguide at E > Umax (case b). 

with E < U m ax gets reflected back; so, these values of E, 
N belong to the classically allowed region. In the oppo- 
site case E > U ma x the particle overcomes the effective 
potential which means that the reflection process is classi- 
cally forbidden. Thus, the line E = Ny/l + a^e" 1 is the 
boundary of the classically allowed region at E, N -C 1. 

When the values of the parameters E, N approach 
the boundary of the classically allowed region the par- 
ticle spends more and more time around the tops of 
the effective barrier U(x). For the values precisely at 
this boundary there are two unstable solutions x — ±1. 
In the two-dimensional picture these solutions corre- 
spond to periodic oscillations around the fixed points 



{x = ±l,y 



a e 



-1/2 



) on the axis of the waveguide. 



We will see that such unstable periodic solutions exist 
beyond the adiabatic approximation and play a key role 




FIG. 3: The boundaries Nt(E) of the classically allowed re- 
gions plotted for ao = 0.5; 0.8; 1.0. The classical reflections 
are forbidden at N < N (E). The points to the left of the 
dashed line E = N correspond to the kinematically forbidden 
initial conditions. 



ripe: \' 



311. Bor- 



in the semiclassical analysis, cf. Refs. [25,, _ 
rowing the terminology from gauge theories [32| we call 
them excited sphalerons" 1 . The sphaleron living at x sa 1 



(x rs —1) will be referred to as near (far) sphaleron ac- 
cording to its position relative to the right end of the 
waveguide. 

To identify the boundary of the classically allowed re- 
gion beyond the adiabatic regime one resorts to numeri- 
cal methods. We scan through the range of initial phases 
< 0o < 2-7T at fixed E, N and look for the reflected 
trajectories. If such a trajectory is found at some </>o, the 
point (E, N) is identified as belonging to the classically 
allowed region. Otherwise the point is attributed to the 
classically forbidden domain. The results of these calcu- 
lations are presented in FigO The boundaries Nb(E) of 
the classically allowed region are obtained for three differ- 
ent values of the parameter ao. The classical reflections 
are allowed for the initial data above the boundaries, 
N > N b (E), and forbidden below them, N < N b (E). 
One observes that beyond the adiabatic regime the form 
of the boundary N (E) depends qualitatively on the value 
of ao. At small ao it monotonically increases. At ao ~ 0.5 
a dip in the curve N (E) develops around E w 0.6. As 
ao grows further, this dip becomes lower and more pro- 
nounced, until it touches the line N = at ao ~ 1. At 
even larger values of ao the boundary of the classically 
allowed region splits into two disconnected parts and a 
range of energies around E — 0.6 appears, where the 
classical reflections are allowed even for N = 0. 

Heuristically, one may envision that the form of the 
curve Nb(E) reflects the behavior of the suppression ex- 



4 The word sphaleron is formed from the Greek adjective aipaXepos 
meaning "ready to fall" . 



ponent F(E, N) in the classically forbidden region. In- 
deed, the suppression exponent is zero above the line 
N = Nb(E). As the value of N gets decreased at 
fixed energy E, the function F(E, N) starts growing at 
N = Nb(E). Thus, the deeper the point is in the classi- 
cally forbidden domain, the larger is F, and vice versa. 
According to this reasoning the dip in the curve Nb(E) 
at ao > 0.5 implies that classically forbidden reflection 
are least suppressed at energies E 0.6. Moreover, at 
ao < 1 there is a finite range of occupation numbers, 
< N < Nb{E = 0.6), where the reflection process is 
suppressed at all energies (except for the narrow band 
N < E < N^/l + a^e -1 corresponding to the adiabatic 
regime). For these values of N, the suppression expo- 
nent F considered as function of E is expected 5 to have 
a (local) minimum in the vicinity of E — 0.6. One of the 
purposes of this paper is to test the method of complex 
trajectories in the regime when the minimum of F(E) 
exists; so, we concentrate on the case ao = 0.8. 

Now, we are in a position to investigate the classical 
dynamics of the system (ITT1) in detail. The observations 
we make below are crucial for the subsequent study of 
the classically forbidden reflections. One asks the follow- 
ing question. At given E, N belonging to the classically 
allowed region, there is a non-empty set TZe,n of initial 
phases 4>q, which give rise to classical reflections. What 
is the structure of this set? To answer this question, we 
fix the initial conditions (JT3J) , (fTi|) at ti = and integrate 
the equations of motion until tt = 200 starting from dif- 
ferent values of the initial phase <po ■ In this way, the de- 
pendence of Xf = x(tf) on </>o is obtained, see Fig. 2] The 
negative values of Xf correspond to the classical trans- 
missions through the waveguide, while Xf > represent 
reflections. Thus, 



Ke.n = {<fo\xf(<fa) > 0} 



(18) 



Figure [4] shows that the set 1Ze,n is not connected: the 
intervals of phases corresponding to the reflected trajec- 
tories are intermixed with those representing transmis- 
sions. Moreover, the scaling of the fine structures of the 
function Xf((f>o) reveals self-similar behavior. One con- 
cludes that the set 1Ze,n consists of an infinite number of 
disconnected domains forming a fractal structure. Such 
a complexity is a manifestation of irregular dynamics in- 
herent in our model; this feature is in sharp contrast to 
the situation one observes in completely regular systems 
(see, e.g., the model of Ref. [33l where the analogous set 
consists of a single interval [30(). 

To analyze the nature of irregular dynamics, we con- 
sider the trajectories generated by the initial phases 
which span various connected intervals TZ a - t E N C TZe.n- 



5 We stress that the heuristic arguments about the behavior of 
the suppression exponent will be confirmed by the explicit semi- 
classical and quantum mechanical calculations in the subsequent 
sections. 
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FIG. 4: The dependence of the final particle position on the 
value of the initial phase (j>o for E = 0.6, N = 0.1 (xf axis 
is not to scale). From top to bottom: scaling of the fine 
structures of the function Xf(<f>o) reveals self-similar behavior. 
The corresponding regions in different graphs are marked by 
the boxes of the same color. The trajectories corresponding 
to the marked points are shown in Fig. [5] 



All the classical trajectories from a given interval of 
phases display the same qualitative properties; some fea- 
tures change discontinuously, however, as one goes to an- 
other interval. Let us characterize each trajectory by its 
behavior in the interaction region. To start with, one 
can distinguish a subset of intervals TZj- e,n C He,n 
corresponding to the trajectories which reach the far 
sphaleron, perform several oscillations there, and go out 
of the interaction region back to x — » +oo. We call this 
subset "the main sequence" . The a;(£)-dependence for the 
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first four trajectories from the main sequence is plotted 
in Fig. [5^,, while the corresponding values of the initial 
phase are marked in Fig. [5] by numbers (1 to 4). Let us 



* 4 




-2 1 1 1 1 1 1 

10 20 30 40 50 

t 

(a) 




-2 1 1 1 1 1 1 1 

10 20 30 40 50 60 

(b) 

FIG. 5: The main (a) and secondary (b) sequences of trajec- 
tories for E — 0.6, iV = 0.1. The trajectories correspond to 
initial phases with the same marks as in Fig(4] 

make the important observation: all the plotted solutions 
leave the far sphalcron with approximately the same os- 
cillatory phase. In fact, we have found that this is true for 
any reflected trajectory. Here lies the root of the discon- 
nectedness of 7Ze,n- one cannot transform one reflected 
trajectory into another with a different number of oscil- 
lations at the far sphaleron by continuously changing the 
initial phase. For a solution from the main sequence the 
index j is equal to the number of oscillations at x ~ — 1 . 

The trajectories corresponding to the intervals oHZe,n 
other than those from the main sequence display more 
involved behavior. After oscillating at the far sphaleron 
they move to the near sphaleron, oscillate an integer num- 
ber of periods on top of it, return to the far sphaleron, 
oscillate there once again, etc. As an example we present 
several trajectories from the next-to-the-main sequence 
(with only one return to the far sphaleron) in Fig. [5}d, 
they correspond to the initial phases marked by la to 
3a in Fig. [4] We observed that the motion back and 
forth between the sphalerons can be arbitrarily compli- 
cated giving rise to the aforementioned fractal structure 
of He.n- 

For the semiclassical analysis of the next Section it is 
important to know what happens with the classical re- 
flected solution when the initial phase (f>o approaches the 
end-point 4> en d of an interval lZ a -,E,N C TZe,n- Figured] 



implies that at the ends of the interval the correspond- 
ing trajectory gets stuck in the interaction region. More 
precisely, one observes that as the value of 0o approaches 
4> en d, the classical solution spends more and more time 
at x K 1 before going away to infinity (see Fig. O, so 
that the initial datum (f>o = 4>end gives rise to a trajec- 




-2 1 ' ' 1 

20 40 60 



FIG. 6: The classical trajectory with the value of the initial 
phase 4>o = —0.414 close to the end-point, <j>a — (fiend ~ 10 -3 , 
of the first interval of the main sequence; E = 0.6, N = 0.1. 
Before leaving the interaction region the trajectory gets stuck 
at the near sphaleron for a long time. 

tory which is neither reflected, nor transmitted, but ends 
up at the near sphaleron. Such a trajectory is unstable: 
there exist arbitrarily small perturbations which push it 
out of the interaction region to either end of the waveg- 
uide. As the number of intervals (and hence of their 
end-points) constituting 1Ze,n is infinite, the number of 
the above unstable trajectories is infinite, too. We will 
see in the next Section that the above unstable solutions 
give rise to the tunneling trajectories; thus, the classically 
forbidden reflection process of the next Section provides 
a particular example of chaotic tunneling. 

More insight into the classical dynamics is gained by 
extending the previous analysis to include the variation 
of two initial conditions, (j>o and N, at fixed energy E. 
One is again interested in the structure of the set TZe of 
initial data (eVi , N) which correspond to classical reflec- 
tions. This set is shown in Fig. [7] (E = 0.6). As before, 
it consists of an infinite number of disconnected domains 
TZ a -E- Each of these domains is characterized by the 
way in which the corresponding trajectories travel back 
and forth between the sphalerons. The domains of the 
main sequence with j = 1, 2, 3 are clearly visible in Fig.0 
along with the secondary ones which accompany them. 
Figure [7] enables one to understand what happens with 
the classical reflected trajectories when the occupation 
number N approaches the boundary of the classically al- 
lowed region from above, N — > Nb(E); this corresponds 
to moving to the lower boundary of the set TZe, 

N b (E) = inf N . 

(<t> ,N)eiz E 

One notices that no matter how close N is to this bound- 
ary, the line N = const always has intersections with 
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some domains from the main sequence. In other words, 
the boundary of the classically allowed region is the accu- 
mulation point of the main sequence of domains. On the 
other hand, as N tends to Nb(E), the intersection of the 
line N = const with certain domains oUZe disappears. 
This means disappearance of certain types of trajectories 
at N — > Nb(E). Concentrating on the main sequence, one 
observes that the remaining trajectories are those with 
large indices j. As the latter are equal to the number 
of oscillations at the far sphaleron, one concludes that 
the closer the point (E, N) is to the boundary of the 
classically allowed region, the longer the corresponding 
trajectories get stuck at the far sphaleron 6 . 

We conclude this section with remarks on the implica- 
tions of the above classical picture for the semiclassical 
description of tunneling. It will be shown in Sec. IIVI 
that each isolated domain of the set TZe continues into 
a branch of tunneling trajectories. Thus, the number of 
tunneling solutions at given E, N < Nb(E) is infinite; 
each of these solutions is associated with a certain do- 
main TZ a -E- In particular, the fractal structure oiH-E is 
inherited by the collection of tunneling paths. The above 
property is very different from that of the models with 
completely regular classical dynamics [30| where the tun- 



This property should not be confused with the property de- 
scribed in the previous paragraph, where it was pointed out that 
the trajectories get stuck at the near sphaleron when the initial 
data approach the boundary of a single domain 7?. a; e of H-e- 



neling solution is unique. On the other hand, an infinite 
number of tunneling trajectories seems to be generically 
inherent in chaotic systems [20(. According to the gen- 
eral rules, the suppression exponent F in Eq. (JTJ) is equal 
to the lower bound of the suppressions calculated on all 
the tunneling solutions, 

F(E,N) =MF a (E,N) , (19) 

a 

where the index a marks the solutions. At first glance 
it is not clear how the lower bound (|T9]) can be found: 
one is unable to compute the suppression exponents of 
an infinite set of trajectories. Two observations which 
greatly simplify the problem are as follows. First, it will 
be shown below that the nearer the domain TZ a . e is to the 
boundary of the classically allowed region, the smaller is 
the suppression exponent of the tunneling solution asso- 
ciated with it. Second, the domains of the main sequence 
accumulate near the boundary of the classically allowed 
region. This implies that the suppression exponents Fj 
of the tunneling trajectories from the main sequence 7 de- 
crease with the index j. Moreover, for any tunneling path 
there exist tunneling solutions from the main sequence 
with the same E, N and smaller value of the suppression 
exponent. Thus, in order to calculate the lower bound 
(fTTj|) it is sufficient to consider the trajectories from the 



7 These are the ones originated from the main sequence TZj . g of 
domains. 
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main sequence, compute their suppression exponents Fj , 
and take the limit, 

F(E,N) = lim Fj (E, N) . (20) 

j—KX 

This tactics is implemented in the next section. 



IV. SEMICLASSICAL STUDY 

We now formulate the boundary value problem for the 
tunneling trajectories. The probability 1(5)) of over-barrier 
reflection is described by the complexified solution to the 
classical equations of motion 



SS 



ss 



5x(t) 6y(t) 







(21a) 



obeying the following boundary conditions, 

• 2 • 2 2 

Y+N = E , y- + L-=N, at * -> -oo , (21b) 



at t -> +oo . (21c) 



Imi , limy — > 0, 



Here E, N stand for the (rescaled, see Eqs. ([13])) energy 
and initial occupation number. 8 The initial conditions 
(|21bj) at t — ► -oo can be cast into the form (|12b[) . (fI3)). 
where the initial phase (f>o is now allowed to take complex 
values. Similarly, x(ti) is also complex, in general. The 
boundary conditions (|2 lb[) . (|21c[) have clear physical in- 
terpretation. Equations (|21b[) fix the quantum numbers 
£, J\f of the initial quantum state. Due to the uncertainty 
principle this makes the two conjugate coordinates, 4>o 
and x(ti), maximally indeterminate. Accordingly, in the 
semiclassical picture they become complex-valued. On 
the other hand, the quantum numbers of the final states 
are not fixed, and Eqs. (|21c)l imply that the particle 
comes out in a classical state with real coordinates and 
momenta. 

Note that, generically, the tunneling solution is defined 
along a certain contour in the complex time plane. In our 
case, however, the contour is trivial: it runs along the real 
time-axis. 

It is useful to parametrize the imaginary parts of x(U) 
and 4>q as follows, 



2Imx(ti) = —x(ti)T 
2Im0 o = -T-6 , 



(22) 
(23) 



where T and 9 are real parameters. Then the suppression 
exponent of a given complex trajectory is 



F a = 2ImS Q - ET a - N9 a 



(24) 



8 Note that the boundary value problem 12111 is invariant with 
respect to time translations: if x(t), y(t) is a solution to Eqs. 
12U . then x(t + r), y(t + r) is also solution for any r £ R. We 
fix this ambiguity by requiring Rex(ti) = 10 (cf. Eq. Ill2al l), 



where the two last terms result from the non-trivial ini- 
tial state of the process, while 



Sq 



dt {-xx - yy - (y - a(x)) 2 ) (25) 



is the classical action integrated by parts. The subscript 
a is introduced to remind that there may exist several 
tunneling solutions with given E, N. 

We do not present the derivation of the above bound- 
ary value problem. The logic is completely analogous to 
that of Refs. [3(1 [H, [H, [35|, and an interested reader is 
referred to these papers. The field theory analog of the 
problem (j2"T|) was first introduced in Ref. 0]. 

It is important to remark that the problem (f2"Tj) does 
not guarantee per se that its solutions describe reflec- 
tions. To ensure that this is the case, one supplements 
Eqs. (|2"l"T) with the condition 



Rex 



at t 



-oo 



(26) 



Below we find solutions which satisfy this requirement by 
using the e-regularization method of Refs. [301 ] . 

Let us explain the physical meaning of the initial-state 
parameters T, 6. One can prove the relations (see, e.g., 
Refs. [Hill) 



T : 



d_ 
~dE 



F(E,N) , 



d 



F(E,N), (27) 



which imply that T and 9 are the derivatives of the sup- 
pression exponent with respect to energy and initial os- 
cillator excitation number. One notices that T, 9 can 
be used instead of E, N to parametrize the tunneling 
paths. Then, the solutions with T = correspond to 
the extrema of the suppression exponent with respect 
to energy. These solutions are called "real-time instan- 
tons" , they can be found directly using the method of 
Ref. [T(|. Calculating the value of the functional (f2"4"]l 
on them, one obtains the extremal (notably, minimal) 
values of the suppression exponent at N = const. The 
method of real-time instantons is important in field the- 
ory [To| . where it enables one to calculate the minimal 
suppression of the collision-induced tunneling. One of 
the purposes of this paper is to check the above method 
by the explicit comparison with the exact quantum me- 
chanical results. Accordingly, we pay specific attention 
to the region around the minimum of the suppression 
exponent F(E). 

We solve the boundary value problem (|2~T|) numeri- 
cally with the deformation procedure. Namely, the solu- 
tion with energy E + AE and oscillator excitation num- 
ber N + AN is found by the iterative Newton-Raphson 
method [36| starting from the solution at E, N, which 
serves as the zeroth-order approximation. In this way, an 
entire branch of tunneling trajectories can be obtained 
starting from a single solution and walking in small steps 
in E, N. The details of our numerical technique can be 
found in Refs. (33|, see Refs. [1, Q for the applications in 
field theory. 
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The only non-trivial task of the above approach is to 
find the input for the first deformation. The idea we 
put forward in this paper is to start the procedure from 
the classically allowed region and arrive at the tunnel- 
ing solutions by changing the values of E, N in small 
steps. [Note that the approach of obtaining complex tra- 
jectories from the real ones was also implemented in a 
different context in Ref.[37].] One begins with the ob- 
servation that the real classical trajectories satisfy 9 Eqs. 
(f2Tj) . Naively, one takes a trajectory with E, N from 
the classically allowed region, and decreases the value of 
N with the hope to get the correct tunneling solution 
at N < Nb(E). However, serious obstacles arise on the 
way. First, the classical solution is not unique at given E, 
N; the degeneracy is parametrized by the initial phase 
4>o € TZe.n- This makes the numerical implementation 
of the deformation procedure problematic. Second, sup- 
pose one takes the initial data (N, </>o) belonging to some 
connected domain !Z a \E cIZe and decreases N. It was 
observed in Sec. IIIII that, as the value of N approaches 
the boundary TV = Nb, a {E) of TZ a -E from above, the re- 
flected classical solution spends more and more time in 
the interaction region. At N = Nb tCl (E) it gets stuck at 
x ps 1 forever, merging at this point with the classical 
solutions, which represent the transmissions of the par- 
ticle through the waveguide at N < Nb, a (E). Here lies 
the worst obstruction: the transmitted classical trajec- 
tories obviously solve the boundary value problem (p?Tj) . 
so that the deformation procedure which starts from the 
classically allowed region will produce them at N < Nb, a , 
instead of the correct complex reflected solutions. One 
concludes that the condition ([26]) should be somehow in- 
corporated explicitly into the boundary value problem. 

A method which automatically fixes the asymptotic 
(I2rj)) of the solution is proposed in Refs. [3(J 01 ■ It 
is called e-regularization. Here we briefly describe this 
method concentrating on its application to the problem 
at hand. A more detailed description of the technique 
can be found in Refs. [3(1 HU . One replaces the action 
of the system in Eqs. (f2Tj) , ([24]) with the modified action 

S £ [x]=S[x]+iel\ nt [x}. (28) 

Here e is small and positive, while the functional Ti nt 
measures the time the particle spends in the interaction 
region. The simplest choice is 

T mt = J dtf(x(t),y(t)) , (29) 

where the function / is real and positive at x, y G M, 
and is localized in the interaction region. Otherwise the 
choice of / is arbitrary: the final result is recovered in 
the limit e — > +0 and does not depend on the particular 
form of this function. We use 

/ = 0(1 - x)9(x - 1) , (30) 



Their suppression is obviously zero. 



where 

- 1+c -2,-,. ( 31 ) 

is the smeared ^-function. Note that / peaks at the near 
sphaleron; this choice will be explained shortly. 

One notices two important changes that the substitu- 
tion l]28p brings into the boundary value problem (|21j) . 
First, the degeneracy of the classical solutions is removed 
at e =/= 0. Indeed, any unperturbed classical trajectory 
extremizes the original action functional S[x], so that at 
fixed E, N and e = there exists a valley of extrema 
parametrized by (f>Q. The functional Ti nt [x], however, 
discriminates among all these extrema and lifts the val- 
ley. Correspondingly, at small e > the extrema of the 
regularized action S e [x] are close to the classical reflected 
solutions with dT int /d4>Q = (see Refs. 30] for the de- 
tailed discussion). From the physical viewpoint this can 
be understood as follows. In the e-regularized case the 
suppression, Eq. of the trajectories in the classi- 

cally allowed region is not precisely zero because of the 
complex term ieTi nt in Eq. (|28[) . The suppression is min- 
imized on the real classical trajectories corresponding to 
the minima of Ti nt \ therefore, the solutions to Eqs. (|2ip . 
(j28|) . which by construction correspond to the least sup- 
pressed reflections, are close to the above classical trajec- 
tories. 

In practice, one starts with the function Ti nt (4>o) rep- 
resenting the values of Ti nt on the classical trajectories 
with fixed E, N and e = 0, finds the extrema of this 
function, and uses the corresponding trajectories as the 
zeroth-order approximation to the solutions of Eqs. (|21[) . 
(f2"8")l with small, but non-zero e. 

The second useful property of the e-regularization 
comes about as one tries to decrease the value of the 
initial oscillator excitation number, and cross the bound- 
ary N — Nb, a (E) of a given classically allowed domain. 
One discovers that at e ^ each reflected trajectory at 
N > Nb, a (E) is smoothly connected with the complex re- 
flected solution at N < Nb, a {E). The reason, again, lies 
in the additional suppression caused by the term ieTi nt in 
the action. Although the reflected trajectories get only 
slightly perturbed at small e, the solutions ending up 
in the interaction region at t — > +oo change drastically. 
Indeed, the functional T;„ t [x] diverges on the latter tra- 
jectories giving rise to the infinite suppression. The im- 
mediate consequence is that the paths which get stuck in 
the interaction region are excluded from the set of solu- 
tions to the boundary value problem (|2"Tj) , (f2"8"|) . Now, the 
solutions to Eqs. (|2lj) . (|28| cannot change their asymp- 
totic (|2l)l): starting from the classical reflected solution 
and decreasing the value of the initial excitation number 
N in small steps, one leaves the classically allowed region 
and obtains the correct complex reflected trajectories. 

Let us briefly comment on the choice of the contour in 
the complex time plane which carries the tunneling solu- 
tions. Obviously, the classical trajectories, which are the 
starting point of our deformation procedure, run along 
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the real time axis. Generically, while moving into the 
classically forbidden region, one may need to deform the 
time contour in order to avoid the singularities of the so- 
lution, see Refs. [3(| ■ In the actual calculations of this 
paper we did not encounter such a necessity. The time 
contour was always kept coincident with the real axis. 

An example of a regularized tunneling solution deep 
inside the classically forbidden region (E = 0.6, N = 0) 
is shown in Fig.[8j3. The depicted solution descends from 
the classical reflected trajectory with E — 0.6, N = 0.1, 
cf>o = —0.366 (see Fig. [5Ji), which (locally) minimizes 
the function Ti nt {4>o)] the regularization e = 10~ 6 was 
switched on and N was gradually decreased by small 
steps. Note that the solution in Fig. [5b is genuinely com- 
plex and satisfies the requirement (|26p. 

The solutions to the original boundary value problem 
are recovered in the limit e — > +0; the suppression expo- 
nents of the unperturbed trajectories are 



F a {E,N)= lim F a e (E, N) 



(32) 



Practically, the removal of the regularization is carried 
out by taking e to be small enough, e < 10~ 6 . At these e, 
the values of the suppression exponents F a ^ e stabilize at 
the level of accuracy O(10 -5 ) which is sufficient for our 
purposes. 

As one tries to take the limit e — > +0 of the tunnel- 
ing trajectory itself, a surprise comes about. Namely, as 
e gets decreased, the regularized tunneling trajectories 
deep inside the classically forbidden region stay longer 
and longer at finite x, so that the trajectories with e = 
do not escape to infinity, but end up on an unstable solu- 
tion - the near sphaleron - living at x sa 1. [An example 
of unperturbed tunneling trajectory (E = 0.6, N = 0, 
e = 0) is shown in Fig. [5b.] Strictly speaking, the so- 
lutions at e = do not describe reflection, but rather 
creation of the unstable state, the sphaleron. Neverthe- 
less, their suppression exponents are relevant for tunnel- 
ing as the latter state decays 10 into the asymptotic region 
x — ► +oo without exponential suppression. The mecha- 
nism of dynamical tunneling via creation of unstable pe- 
riodic orbits was discovered independently in Refs. [301 ] 
and Refs. [H|. 

We remark that all unregularized tunneling trajecto- 
ries in the model under consideration tend to the unstable 
sphaleron solution at t — > +oo, and thus turn out to be 
unstable themselves. Straightforward numerical methods 
are inappropriate for finding such trajectories [331] . while 
e-regularization provides a universal method for treating 
such instabilities 11 . 



10 Classically, the particle sits at the sphaleron for an infinite period 
of time. However, quantum fluctuations lead to the decay of this 
state with the carachteristic time of order \ng. 

11 An alternative would be to change the boundary conditions 
H21cl l, see App.fAl 




t 

(b) 




FIG. 8: (a) The reflected classical trajectory from the first 
domain of the main sequence, corresponding to the minimum 
of T int {4> Q ) (E = 0.6; N = 0.1; <j> = -0.366). (b) The 
regularized tunneling solution (E — 0.6, N — 0, e — 10~ 6 ) 
descended from the above classical trajectory, (c) The limit 
e — > of the tunneling solution. 



Let us now discuss the structure of tunneling solutions. 
Following the above strategy, one starts with the classi- 
cal reflected solutions. It is straightforward to see that 
the function Ti nt (4>o) has at least one extremum in every 
connected interval TZ a - e,n in the set of reflection phases 
1Ze,n- Indeed, the classical reflected trajectories with 
4>a G 7la;E.N spend finite time in the interaction region 
and smoothly depend on the initial data, producing the 



smooth function Ti nt (4>o). On the other hand, the classi- 
cal trajectories spend more and more time at the near 
sphaleron as <j>o approaches the end-points of lZ a ;E,N 
(see Sec. Thus, Tj nt ((/>o) tends to infinity at the 

end points of this interval, necessarily attaining a min- 
imum somewhere in between 12 . [Note that the specific 
form of the functional Ti nt (function / in Eq. (f2"9"| peaks 
at a; i=a 1) has been chosen in accordance with the ten- 
dency of classical trajectories to get stuck at the near 
sphaleron.] Thus, starting from each interval TZ a ;E,N, 
one obtains one branch of solutions to the e-regularized 
problem. Note that, as TZe.n is the section of the set 
TZe by the line ./V = const, the domains of TZe may cor- 
respond to two disconnected intervals of IZe.n, see Fig. [7] 
When TV gets decreased, these intervals merge. Below we 
assume that the starting classical solution of the type a 
has been taken at small enough TV, such that the domain 
TZ a -,E gives rise to a single interval of phases 13 TZ a -E,N- 
Then, the fractal structure of TZe is inherited by the com- 
plex tunneling paths: the distinct branches of complex 
trajectories are in one-to-one correspondence with the 
connected domains of the set TZe- 

Now, we can give a heuristic justification of the claim 
made in the end of Sec. IIIII the nearer the domain 
TZ a; E C TZe lies to the boundary of the classically for- 
bidden region, the less suppressed is the corresponding 
tunneling trajectory. Let us consider two branches of so- 
lutions stemming from the domains 1 and 2 in Fig. [7] 
The line TV = 0.07 intersects with the domain 2, so that 
the point E = 0.6, TV = 0.07 belongs to the classically 
allowed region, and the suppression of the solution 2 van- 
ishes in the limit e — » 0. On the other hand, the line 
TV = 0.07 does not cross the domain 1. Thus, at E = 0.6, 
TV = 0.07 the corresponding solution to Eqs. $ZT$, fl25J| is 
genuinely complex and has non-zero suppression even at 
e — > 0. Suppose that one decreases TV and enters the clas- 
sically forbidden region. At some point the solution of 
the type 2 also becomes classically forbidden. It is clear 
that at least for some range of TV, its suppression remains 
weaker than that of the solution 1. This suggests, though 
does not guarantee, the same hierarchy of suppressions 
inside the entire classically forbidden region. Below, we 
check the conjectured hierarchy explicitly. 

Let us discuss in detail the case TV = (no transverse 
oscillations in the initial state) corresponding to the ex- 



In general, there might be several extrema of Ti n t inside each 
connected interval "R- a -E N> an< i one should consider solutions 
corresponding to each of these extrema. In our case the minimum 
is unique. 

One legitimately asks what happens if the starting classical solu- 
tion is taken at large enough N, where two different intervals of 
phases correspond to one and the same domain of TZe ■ Clearly, in 
this case the above deformation procedure will produce two dif- 
ferent solutions to Eqs. ( 1211 . (I2S D . We have observed, however, 
that as the value of N gets decreased below the point where the 
two intervals merge, the corresponding solutions become identi- 
cal. 
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1 


0.3188 


0.1625 


0.1098 


0.1398 


0.2259 


2 


0.2586 


0.1380 


0.0991 


0.1341 


0.2221 


3 


0.2373 


0.1340 


0.0979 


0.1336 


0.2219 


4 


0.2307 


0.1333 


0.0978 


0.1336 


0.2219 


5 


0.2285 


0.1331 


0.0978 


0.1336 


0.2219 


oo 


0.2272 


0.1331 


0.0978 


0.1336 


0.2219 



TABLE I: The suppression exponents of the complex trajec- 
tories from the main sequence at TV = 0. The rows represent 
the indices j of the trajectories, while the columns correspond 
to the values of energy E. The last row refers to the limiting 
solution. 

treme values of parameters inside the forbidden region; 
all the qualitative features are the same for other val- 
ues of TV as well. Our strategy is to find the limit (|20|) 
of the main sequence of suppressions Fj, and check that 
this limit represents the lower bound of the suppression 
exponents of all the tunneling solutions. 

Figure [9] shows several first representatives of the main 
sequence of tunneling paths at E = 0.5, TV = 0, e = 10~ 6 . 
Similar to their classical progenitors, the branches of 
tunneling trajectories may be classified according to the 
number of oscillations they perform at the far sphaleron. 
On the contrary, the number of oscillations at the near 
sphaleron is not an invariant of the branch: as discussed 
above, it grows as the rcgularization parameter e gets 
decreased, so that the tunneling trajectories get stuck at 
the near sphaleron at e = 0. Another important observa- 
tion about the plots in Fig. [5] is that the imaginary parts 
of the solutions are sizeable only in the beginning of the 
evolution. They fall off rapidly during the oscillations at 
the far sphaleron and become small at late times. This 
qualitative feature holds for the trajectories at other en- 
ergies as well 14 . Thus, the suppression exponents of the 
trajectories are saturated during the first few oscillations 
at the far sphaleron, and depend weakly on the subse- 
quent evolution. In addition, the trajectories from the 
main sequence almost coincide in the beginning of the 
process. This implies fast convergence of the suppression 
exponents Fj to the limiting value F according to the 
formula (|2"0"]). 

The above convergence is demonstrated in Table HI 
where the suppression exponents Fj (E) are presented for 
several values of energy E at TV = 0. As expected from 
the heuristic argument, the value of Fj decreases as j 
gets larger. The limiting values F(E) are also included 
into Table [H They can be obtained by extrapolating the 
dependences of the suppression exponents Fj on j, which 



On the contrary, the smallness of Imi at t —* — oo observed in 
Fig. [9] is peculiar to the trajectories at E = 0.5. It implies that 
the corresponding values of T are small, see Eq. [1 2 2 [ > . so the 
plotted trajectories are close to the real-time instantons. 
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FIG. 9: Four tunneling solutions from the main sequence; 
E = 0.5, N = 0, e = 1(T 6 . 




FIG. 10: The limit of the main sequence of tunneling tra- 
jectories: the solution which gets stuck at the far sphaleron; 
E = 0.5, N = 0. 



are well fitted by the formula 

Fj= F + ae- bj , 

where a, b are real positive coefficients. A better way, 
which is exploited in this paper, is to find the limit j — > 
+oo of the tunneling solution itself, and then calculate 
the limiting value of the suppression exponent using Eq. 
(|24[) . The limiting solution performs an infinite number 
of oscillations on the far sphaleron; it can be obtained 
numerically by the method described in App. A. The 
limiting solution at E = 0.5, N = is shown in Fig. [TU1 

So far we have considered only the tunneling trajec- 
tories from the main sequence and demonstrated that 
Eq. ([20")) reproduces the lower bound of their suppres- 
sions. We have checked that the limiting value F(E) is 
lower than the suppression exponents of other tunnel- 
ing trajectories as well. As an illustration, consider the 
tunneling trajectories shown in Fig. 111! The differences 
between their suppression exponents and that of the lim- 
iting solution are plotted in Fig. 1121 The limiting solution 
is evidently the least suppressed one. 

Let us discuss the physical interpretation of the ob- 
tained results. The form of the limiting solution im- 
plies that the reflection process in our model proceeds 
in two stages. First, the far sphaleron state gets cre- 
ated. This stage is exponentially suppressed due to the 
essential modification of the particle state needed for the 
sphaleron creation. Second, the sphaleron state decays 
into the asymptotic region x — > +oo with probability 
of order 1. This two-stage process is a manifestation 
of the phenomenon of "tunneling on top of the barrier" 
[30L |3"H ; the phenomenon is generic for the inclusive tun- 
neling processes, see Refs. [3, H3] for the field theoretical 
examples. 

Remarkably, though we came to the limiting solution 
considering the accumulation point of an infinite num- 
ber of complicated tunneling paths, the solution itself is 
unique and very simple. All the chaotic features of mo- 
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FIG. 12: The difference between the suppression exponents of 
the secondary trajectories of the types shown in Fig. [11] and 
that of the limiting solution; N = 0. 



tion, like going back and forth between the sphalerons, 
are related to the second stage of the reflection process: 
the decay of the far sphaleron. As the second stage is un- 
suppressed, one might conclude that, after all, the chaotic 
motions are irrelevant for the calculation of the main sup- 
pression exponent. Note, however, that one is not able to 
guess from the beginning, which tunneling solution corre- 
sponds to the smallest suppression; therefore, the system- 
atic analysis of the chaotic motions is essential. Besides, 
the sub-dominant contributions of the chaotic trajecto- 
ries might be important for the analysis of the process at 
finite values of the semiclassical parameter g [l3|, EH HH • 
The final semiclassical results for the dependence of 
the suppression exponent F on energy E for N = 
0; 0.02; 0.04 are presented in Fig.UHl Note that the sup- 
pression exponent is non-monotonic function of energy 
with the minimum near E — 0.5. 15 We remind that the 
minima of F(E) correspond to the particularly interest- 
ing solutions, real-time instantons. 



V. EXACT QUANTUM COMPUTATIONS 

In this Section we extract the suppression exponent F 
from the exact reflection probability ©. One begins by 
solving numerically the Schrodinger equation 



H|V)=£|V) 



(33) 



It is convenient to work in the original variables X, Y , 
Px , Py , see Eq. (J5j) , in order to bring the kinetic term of 
the Hamiltonian into the canonical form. One rewrites 



FIG. 11: The examples of secondary tunneling trajectories; 
E = 0.5, N = 0, e = 1(T 6 . 



15 In the range E < 0.1, which is not shown in Fig. 1131 the sup- 
pression attains the local maximum and goes down to zero as E 
decreases toward the boundary of the classically allowed region. 
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FIG. 13: The semiclassical (lines) versus quantum mechanical (points) results for the suppression exponent F(E) at three 
different values of N. The errors of the quantum mechanical computations are smaller than the point size for E > 0.3. 



Eqs. ©, © as 

p2 , p2 i 

H = ix+iV + l (y _ a(x))2) (34) 

where 

a(X) = -a(gX) . (35) 
9 

Basically, our numerical technique follows the lines of 
Refs. [33,[3^|. One works in the asymptotic basis formed 
by a direct product of the translatory coordinate eigen- 
functions \X) and eigenfunctions \n) of the y-oscillator 
with the fixed frequency oj = 1, 

i> n {X) = (X, n\il>) • 
The stationary Schrodinger equation (|33[) reads, 

^MX) =J2 A nn'(X)^(X) , (36) 

where 

A nn ,(X) = (nlP^ + (y -~a{X)f ~2£\n') 

= 8 n ^{2n+l + a 1 {X) - 2£) 

- a(X)(8 n , n > -i\/2n + 2 + 5 n ^ + iV2n) 

is an infinite three-diagonal matrix. 

In the asymptotic regions X — * ±oo, a(X) — > the in- 
teraction terms become negligibly small, and the solution 
takes the form, 

MX)~^r^ x +tte-^ x , (37) 



where 

P n = v/2£ - (2n+T) (38) 

stands for the asymptotic translatory momentum of the 
n-th mode. The boundary conditions for the station- 
ary wave function ip n (X) are constructed in the standard 
way. The particle comes from the right, X — » +oo, in the 
Af-th oscillator state; hence, we fix 

t+ = SnM ■ (39) 

On the other hand, only the outgoing wave should remain 
at X — > — oo, 

r- = . (40) 

Note that while the low-lying modes are oscillatory at 
the asymptotic, the ones with n > £ — 1/2 grow (de- 
cay) exponentially, see Eq. {35]). Physically, the latter 
correspond to the kinetically inaccessible region, where 
the energy of the transverse oscillations exceeds £. We 
fix the boundary conditions for them by killing the parts 
growing exponentially toward infinities. One notes that 
after the proper continuation of Eq. (|38[) , 

P n = i^/(2n + 1) - 2£ , n>£-l/2, (41) 
the aforementioned conditions coincide with Eqs. (|39[) . 

go]). 

Equations (|36|) . (|39|) . (|40|) constitute the boundary 
value problem to be solved numerically. After the sta- 
tionary wave function ip n {X) is found, one calculates the 
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probability current 



J = Im^<(X)^Vn(AM 



and hence the reflection probability 

|j(out)| ^ p n 



V 



IJ( 



in) I 



E 



~+|2 



n<£-l/2 



(42) 



(43) 



where the currents and j( out ) are computed with 

the incoming and outgoing parts of the wave function at 
X — > +oo, respectively. 

The details of the numerical formulation of the prob- 
lem pp]). ([55)1 , [gDD are presented in App. [B] Let us 
discuss the results. The typical dependence of the re- 
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FIG. 14: The function V(E) plotted for N = 0, g = 0.2. 
dashed curve represents the smoothened result, Eq. (|45 



The 



flection probability "P on energy E — g 2 £ is shown in 
Fig. [TU solid line (the values of the other parameters are 
N = g 2 M = 0, g = 0.2). The striking feature is that the 
graph is modulated by sharp oscillations. At first glance, 
this picture is incompatible with the results of the semi- 
classical analysis, where we obtained, at least in the lead- 
ing order, that the tunneling probability, V cx e~ F ^ E ^ 9 , 
is a smooth function of energy. We are going to demon- 
strate the opposite: the quantum mechanical results rec- 
oncile nicely with the semiclassical ones. 

By computing the reflection probability at different g, 
one notices the following important properties. First, the 
period of oscillations scales like te cx g 2 . So, in the semi- 
classical limit g — > the oscillations become more and 
more frequent. Second, one considers F g = g 2 h\V and 
asks whether the amplitude of oscillations of this quantity 
goes to zero in the limit g — > 0. We certainly observed 
that it drops down as g decreases though we were unable 
to figure out 16 whether it indeed vanishes at g — ► 0. The 



properties mentioned above suggest the interpretation of 
the oscillations as the result of the quantum interference 
between the contributions of various tunneling trajecto- 
ries found in Sec. HVl Suppose for simplicity that there 
are only two such trajectories. Let us denote their com- 
plex actions by S\{E) and S 2 (E) (we omit the depen- 
dence on iV). Then, the total reflection amplitude reads, 



*Si(£)/9 2 



+ A 2 e 



iS 2 (E)/g 2 



For the sake of argument we suppress the index n and 
disregard the initial-state contributions into the expo- 
nents. In the above formula A\ and A 2 stand for the 
pre-exponential factors of the partial processes. For the 
reflection probability one writes, 



-21m S 2 (E)/g 2 



V cx |rf - l^xfe 

+ 2\A 1 A 2 \e~ I ™( s ^ +s ^ E »/° 2 x 
x cos [Re (Si(E) - S 2 (E))/g 2 + arg^/^)] . (44) 

Along with the terms corresponding to the probabilities 
of the partial processes, one gets the interference term, 
which results in oscillations with period te of order g 2 . 
Of course, if one of the solutions, say, the first one, is 
dominant, ImSi < ImS^, the relative contributions of 
the two last terms in Eq. (|4~4"|) vanish exponentially fast 
at g — > 0. In our case the situation is more subtle, how- 
ever. As the semiclassical analysis reveals, in the model 
under consideration there exists an infinite number of 
tunneling paths, which pile up near the dominant limit- 
ing solution. So, at each finite value of g there are solu- 
tions which satisfy Im (52 — Si) < g 2 (the index "1" still 
marks the dominant solution here). Thus, the sum (14~4"|) 
always contains an infinite number of oscillating terms 
producing a complicated interference pattern; it is not 
clear whether the oscillations disappear at small g. 

What saves the day is the aforementioned scaling of the 
oscillation period. Indeed, te vanishes at g — * implying 
that the oscillations become indiscernible in the semiclas- 
sical limit. One obtains a quantity which is well-behaved 
in this limit by averaging the reflection probability over 
several periods of oscillations. To be more precise, we 
consider the smoothened probability 



V (s) {E,N) = J dE'V a (E - E')V(E',N) , 



(45) 



where T> a is the bell-shaped function, 



D a {E) 



dED a (E) = 1 



If a = <? 2 £, where S is a fixed number, the smoothening 
(|45|) does not spoil the value of the dominant suppression 



This is due to the limitations of the numerical approach: one 



cannot obtain solutions to the Schrodinger equation at arbitrarily 
small g, see App. [B] 
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exponent. Indeed, for the first term in Eq. (|4"4"|) one 

writes, 

V (s ' dom) cx J dE'lArfe-^'Vo'v^E' - E) (46) 

« (Aipe^ 2 ^^ 2 / 4 • e" Fl(B)/s2 , 

where 17 -Fi = 2ImSi, = dFi/dE, and we made use 
of the Taylor series expansion in the second equality. It 
is clear that only the pre-exponential factor is affected 
by the integration (|45|) , while the exponent Fi (E) is left 
intact. On the other hand, at large enough values of the 
coefficient £ the formula (|45[) represents averaging over 
many oscillatory periods, which kills all the oscillating 
contributions. Indeed, consider the typical interference 
term, 

■P (osc) cx C{E)e' Fo{E]/9 \ 

x cos [AS(E)/g 2 - A<P(E)] , (47) 

where F = Im(Si + S 2 ), AS = Re(S*i - S 2 ), C = 
2\AiA2\, A0 = arg(^4i/A2). Performing integration, one 
obtains, 

p(«,«c)(j0) oc Ce- Fo/92 cos [AS/g 2 +Ac/>- i^AS"S 2 /2] 
xexpj^F^-AS' 2 )} . 

Taking into account that Fq « F[ we see that the in- 
terference term is suppressed by the additional factor 
exp{— £ 2 AS" 2 /4} with respect to the dominant contri- 
bution (|46|) . Below we fix S = 1. We observed that 
the interference patterns get multiplied by 10 -3 in this 
case; the latter number is accepted as the precision of 
the smoothening. The graph of the function p( s \E) is 
shown in Fig. [HI dashed line. 

Our final remark concerns the physical meaning of the 
smoothening procedure. In a realistic experiment one 
cannot fix the energy £ of the incoming particles exactly; 
rather, one works with some sharply peaked energy dis- 
tribution T> of a width A£ ~ E. Formula (|4"5|) represents 
averaging over such distribution. 

Now, we are ready to consider the limit g 2 — > 0. Our 
aim is to check the following asymptotic formula, cf. 
Eq. Q, 

V {s \E,N) -> g^A{E,N)e- Fi - E ' N)/92 as g 2 -> . (48) 

We compute the value of the quantity — g 2 \nV^ at sev- 
eral j 2 « 1, keeping E and N fixed, and fit the graph 
with the expression 

- 5 2 ln-p (s) = F-~/g 2 \ng-g 2 \nA. (49) 



For the sake of argument we, again, disregard the boundary terms 
in the suppression exponent, c.f. Eqs. 11241 . 
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FIG. 15: The quantity -g 2 lnV (s) viewed as a function of q 
at fixed N and E = 0.5. 

It is important to point out that we do not consider 7 as 
a free parameter of the fit. Rather, we use the following 
result of Ref. (34|: for the sphaleron-mediated processes, 
such as ours, 7 = 1 at N — (vacuum initial state), and 
7 = 2 at N 7^ 0. Some of the numerical results (points) 
together with their fits by the formula ([49]) (lines) are 
shown in Fig. 1151 The graphs are drawn for three dif- 
ferent values of N and E = 0.5. Apparently, our data 
are well approximated by the asymptotic (|49|) . To get 
the value of the suppression exponent F, one extrapo- 
lates the curves in Fig. |T5] to the point g 2 = 0, where 
the quantum mechanical results should coincide with the 
semiclassical ones. The error of the extrapolation arises 
mainly from the disregarded terms proportional to the 
higher powers of g 2 . One reduces such errors by comput- 
ing the reflection probability at the smallest possible g 2 . 
In practice we used two values of the semiclassical pa- 
rameter in each fit, namely, g — 0.1, 0.07 at N — 0, and 
g = 0.1, 0.08 at N = 0.02, 0.04. The extrapolation error 
is determined by pouring some additional points into the 
fit; it varies from SF ~ 10~ 2 at small E to SF ~ 10~ 4 
at E ~ 1. It matches with the estimate 0{g i ) ~ 10~ 4 of 
higher-order terms of the semiclassical expansion which 
are neglected in Eq. (|4*9")) . 

Our final results for the suppression exponent F(E, N) 
extracted from the numerical solution of the Schrodinger 
equation are presented in Fig. [T^] (points with error bars 
standing for the accuracy of the extrapolation). The 
quantum mechanical results are in very good agreement 
with the semiclassical ones (lines) . This justifies the semi- 
classical approach presented in this paper. 



VI. SUMMARY AND DISCUSSION 

In this paper we tested the semiclassical method of 
complex trajectories in the regime of chaotic dynamical 
tunneling. We studied a particular example, over-barrier 
reflection in the two-dimensional waveguide model 
The initial state of the process was fixed by the total en- 
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ergy £ and occupation number M of the transverse oscil- 
latory motion. We calculated the suppression exponent 
of the process both semiclassically and by solving exactly 
the full Schrodinger equation. The two approaches show 
very good agreement. 

The tunneling trajectories in our semiclassical ap- 
proach are obtained as solutions to the boundary value 
problem (JUJ). We advocated a particular method for 
finding these solutions. It consists of two important in- 
gredients: 

(i) determination of the solution at some values of the 
initial-state parameters E, N; 

(ii) gradual deformation of the solution to other E, N. 
The deformation procedure (ii) was developed in Refs. 

[1, [HI]. The advantages of this procedure are the sim- 
plicity of numerical implementation and generality. It 
can be applied efficiently to systems with many degrees 
of freedom, including non-trivial models of field theory, 
see Ref. @. 

On the other hand, a generic approach for performing 
the step (i) was missing so far. In this paper we pro- 
posed the systematic procedure which fills this gap. The 
procedure enables one to obtain the complex tunneling 
trajectories starting from the real classical solutions. It 
is based on the e-regularization method of Refs. [30l ] . 
Our procedure appears to be generic. It can be applied 
to any process which proceeds classically at some values 
of the initial-state parameters (at large N in our case) 
and becomes exponentially suppressed at other values. 

The process we studied is a particular example of 
chaotic tunneling. The chaoticity manifests itself in the 
infinite number of tunneling solutions. Our procedure, 
which connects the tunneling solutions to the classi- 
cal ones, turns to be highly efficient in this situation. 
Namely, it enables one to classify the tunneling trajecto- 
ries on the basis of the analysis of the classical dynamics. 
This was demonstrated explicitly in the present paper. 

In addition, we proposed a heuristic criterion for sort- 
ing out the least suppressed tunneling trajectories basing 
on the classification of their classical progenitors. Hope- 
fully, this criterion will be useful for the processes in other 
dynamical systems as well. 

Another interesting feature of our setup is the phe- 
nomenon of optimal tunneling. Namely, the suppression 
exponent F considered as a function of energy E at fixed 
N is non-monotonic. It attains the local minimum at 
E ps 0.5 which is thus (locally) the optimal energy for 
tunneling. It is worth stressing that this behavior of the 
suppression exponent is unrelated to the quantum inter- 
ference, which was neglected in the semiclassical anal- 
ysis and eliminated from the exact quantum computa- 
tions. Another example of tunneling process with non- 
monotonic dependence of the suppression exponent on 
energy is considered in Ref. [35[ . 

Let us mention some open issues. 

The method of complex trajectories is known to suffer, 
in general, from the difficulties related to the Stokes phe- 
nomenon [12l . 1 1 31 ] . The essence of this phenomenon is that 



solutions of the tunneling problem may be unphysical in 
some regions of the parameter space; the contributions 
of such solutions to the tunneling probability should be 
dropped out. Presently, there are no generic criteria for 
distinguishing between the physical and unphysical so- 
lutions. One notes that the method we put forward in 
this paper trivially excludes some of the unphysical solu- 
tions. Indeed, we start from the classically allowed region 
of initial data (large N) where the physical solutions are 
precisely the real- valued classical trajectories. At the 
second step we relate these solutions to the complex tun- 
neling trajectories; thus, we exclude the branches of un- 
physical solutions which are complex-valued deep inside 
the classically allowed region of initial data. Due to this 
"automatic" criterion, we did not see any manifestations 
of the Stokes phenomenon in the semiclassical calcula- 
tions presented above. We remark, however, that the 
tunneling process we studied is the simplest one from the 
point of view of the Stokes phenomenon; the above "au- 
tomatic" exclusion of unphysical solutions is insufficient 
in somewhat more involved situations. Two aspects of 
our model simplify the analysis. First, we observed that 
the tunneling solutions obtained from the real classical 
trajectories of a given topology a form a single smooth 
branch which covers the entire range of initial data. Sec- 
ond, we showed that the hierarchy of the suppression ex- 
ponents F a is the same at different values of E, N. These 
two observations allowed us to identify a single smooth 
branch of physical solutions which give the dominant con- 
tribution to the tunneling probability. In other models 
(see, e.g., Ref. [35J) the solutions obtained by small de- 
formations from different parts of the classically allowed 
region of initial data may correspond to different smooth 
branches of complex trajectories. Each of these branches 
may be physical and dominant in some region of the ini- 
tial data plane, and unphysical or sub-dominant in other 
regions. If this is the case, some method of treating the 
Stokes phenomenon [3, 0, EH should be exploited. 

The approach adopted in this paper was to perform 
the evaluation of the tunneling probability as the sys- 
tematic semiclassical expansion in terms of g 2 . We have 
calculated the leading term (suppression exponent). It is 
of interest to develop a method for calculating the sub- 
leading terms, in particular, the pre-exponential factor 
A, see Eq. ( [IT). Presumably, this can be done along the 
lines of Ref. [34]]. An important problem which should 
be solved here is how to cope with the infinite number 
of tunneling paths, whose suppressions can be arbitrar- 
ily close to the limiting value. In particular, one has 
to understand whether the contributions of all the paths 
should be summed up or the correct value of the prefac- 
tor is determined by a sort of limiting procedure similar 
to the one we used to obtain the suppression exponent. 

Another question is the following one. While solving 
the Schrodinger equation we observed that at finite g 2 the 
dependence of the exact tunneling probability on energy 
is modulated by oscillations. We conjectured that they 
result from the quantum interference of different tunnel- 
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ing trajectories. In the true semiclassical limit g 2 — ► the 
above oscillations become infinitely frequent and should 
be averaged over. In real physical situations one deals, 
however, with small but finite values of the semiclassical 
parameter g 2 . Thus it would be interesting to reproduce 
the interference pattern mentioned above in the semi- 
classical approach. As suggested by Refs. 0, H3, HH 
this could be done by summing up the contributions of 
various tunneling trajectories at finite g 2 . We leave this 
investigation for the future work. 
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the following way. One finds a solution to the original 
equations of motion which spends finite time on the far 
sphaleron, and cuts it at the moment t — tf when it is 
still at x w — 1. Using this trajectory as the zeroth-order 
approximation, one applies the Newton-Raphson algo- 
rithm and obtains the solution satisfying the boundary 
conditions (|A3[) . The latter solution is defined inside the 
interval t £ [tj , tf] . The next step is to continue the solu- 
tion to a larger time interval by gradually increasing the 
value of tf and deforming the tunneling solution. As a 
result, one obtains the solution living at the far sphaleron 
for an arbitrarily long time. 

The solution obtained in this way does not, strictly 
speaking, satisfy the boundary conditions (|21cp of the 
tunneling problem. In order to restore the original 
boundary conditions (|21c[) one should, in principle, in- 
vestigate the dependence Xf{Xf) and find the value of 
x® which solves the equations Xf{Xf) — Xf. However, it 
is not necessary in our case: at late times the limiting 
solution is almost real, so one obtains \mxf ~ auto- 
matically. 

APPENDIX B: NUMERICAL SOLUTION OF 
THE SCHRODINGER EQUATION 



APPENDIX A: A METHOD TO OBTAIN THE 
UNSTABLE SOLUTIONS 

In this Appendix we describe the numerical method 
used to obtain the limiting tunneling trajectories consid- 
ered in Sec. IIV1 namely, the ones ending up at the far 
sphaleron at late times. The main problem here is the 
instability of the trajectories in question. 

Our main idea is to add the following term to the action 
functional, 



S[x] h+ S[x] + iM(x{t f ) - x° f ) 



(Al) 



where x(tf) stands for the final value of the coordinate 
x, while M > 0, x® are real parameters. For a given 
trajectory the term (|A1|) leads to additional contribution 
to the suppression exponent, 



AF = 2M(x{t f )- x)y . 



(A2) 



The introduction of the above term induces the following 
modification of the boundary conditions (|21c|) . 



lmx(t f ) = 0, (A3a) 
Imx(i/) = -2M(x{t f ) - x)) . (A3b) 

This can be shown using the systematic approach of Refs. 
[U H, H]. At lar g e positive M the term (fM|) fixes 
the value of the final x-coordinate of the solutions to be 
close to Xf. For our purposes we choose x° f to be in the 



vicinity of the far sphaleron, x° f = — 1. Solutions which 
approach the far sphaleron at late times are obtained in 



Here we present the numerical formulation of the prob- 
lem ([55]). (fMjl . ijiD]). First of all, the range of the transla- 
tory coordinate should be bounded, — L < X < L, as well 
as the oscillator excitation number, n < N y . Besides, we 
introduce the uniform lattice with spacing A, 



.,N X + 1 



X k = -L+(k + l)A, k = -l, 

where N x = — 2 + 2L/A. The Taylor series expansion 
gives, 

1 



A 2 



[ipk+i ~ 2 ^fc + V'fc-i] 



^'+12^' 



(IV) 



0(A 4 ) . (Bl) 



By ip£ and ip k we denote the second and fourth deriva- 
tives of ip at X = X k . Note that the index n is sup- 
pressed: hereafter we use the matrix notations. From Eq. 
([3"6"]l i^k = Ak4>k', for the fourth derivative one writes, 



-^[< +1 -2^' + ^ 1 ]+0(A 2 



= — [A k+1 ip k+1 - 2A k t(; k + A fc _iVfc-i] + 0(A 2 ) 



Substituting the above expressions into Eq. (|B1[) . one 
gets the forth-order Numerov-Cowling approximation 
for Eq. Q5|l, 



A 2 

1 - j2 Ak+1 > ^ k+1 



A 2 



5A 2 

2 + — A k I Vfc + 



1 - -y^fc-i 1 1>k-i = 



fc = 0,...,A^. (B2) 
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It is worth noting that Eq. (|B2p supports the discrete 
probability current, which is conserved exactly, 



j(d) = 1^+ A * Ak ] ( ;, 



12 



This conservation law was used to estimate the round-off 
errors. 

The boundary conditions ([39)) , (|40|) are imposed at the 
very last and first sites, k = N x + 1 and —1 respectively. 
One notes that the asymptotic formula (|3~T)) holds in the 
discrete case as well, provided the continuum dispersion 
relation is replaced with the discrete one 18 , 

Pn -» J P,i d) = |-arcsin — - P " A . 

A 2y/l + P 2 A 2 /12 

Consequently, one rewrites Eqs. ([55]) . (f^D|) as 







(B4) 



These relations together with Eq. (|B2[) form a system 
of N y (N x + 3) linear equations for the same number of 
unknowns ip n ,k- After solving them, one calculates the 
reflection probability by making use of the discrete cur- 
rent, 



V 



|j( d 



,out) I 



IJ( d > in )l 



where j( d < in ) and j( d > out ) are the incoming and outgoing 
currents at k = N x + 1, 



j(d,< 



- y 

n<£-l/2 



sin(P T ( t d) A)|r+| 2 
(l-Isin 2 ^) ; 



j(d,in) = _ 



1 sin(P^ d) A) 



A 



1 _ i S i n 2 £al^ 



2 ' 



while the reflection amplitudes r+ are extracted from the 
wave function, 

/ -iP (d) A / 

" l_ e -2i4 d) A ' 

We will see shortly that the finite difference approx- 
imation works well only if the numbers of the lattice 
points and oscillator levels are large. Typically, N x ~ 
10000, N y ~ 500, and the system pi]), (|B2]) contains 
N y (N x +3) ~ 5- 10 6 equations. Such an enormous system 
of equations cannot be solved with the general algorithms 



For n > £ - 1/2 one has = farcsh -^JM0= = 



of linear algebra. So, we took advantage of the special 
form of Eqs. (|B2[) . Namely, the fc-th matrix equation 
relates the vector ipk to the unknowns at the adjacent 
sites tpk-i and ipk+i only; by performing numerically the 
matrix inversion it can be recast in the form 

ipk = ifcV'fc-i + Rktpk+i , 

where an Rk are the N y x N y matrices. One sub- 
stitutes the above formula into the other equations of 
the system (|B2[) . thus excluding ipk, as well as the k-th 
matrix equation. Performing this operation repeatedly, 
one ends up with a few matrix equations, which can be 
solved in a straightforward manner by the LU decompo- 
sition method. It is worth pointing out that the variables 
ij)k and tp q which are not neighbors to each other, can be 
excluded in parallel, so that the above algorithm is suit- 
able for the multiprocessor machines or computational 
clusters. The reader interested in the details of the al- 
gorithm should address Refs. [Hj], (38[, or our Fortran 
90 code [1^], which hopefully can be executed on other 
machines. 

Before proceeding to the actual calculations, one makes 
sure that the parameters of the lattice are chosen prop- 
erly, so that the truncation and discretization errors are 
kept under control. Our purpose is to get the quantum 
mechanical results in the semiclassical region g 2 — > 0. 
Therefore, it is convenient to account explicitly for the 
dependence of the lattice parameters on g. The trunca- 
tion L of the translatory coordinate is fixed by the con- 
dition that the interaction represented by the function 
a(X) is small enough at X = L. Taking into account the 
scaling (|3"B"|) of a(X) with g, one obtains the formula 

L = L/g, 

where L is fixed and large. In the practical calculations 
we used the value L = 12, which is large enough as 
a(L) ~ 10 -31 : at g > 0.07 this number is smaller than 
the absolute values of the reflection amplitudes, the lat- 
ter exceeding 10 -8 . We have chosen the truncation of the 
oscillator levels in accordance with the condition that the 
occupation number of the last mode is negligible, 



£1* 



W„-l,fc 



< 10 



-30 



this inequality was satisfied with N y typically ranging 
in between 200 and 500. The last parameter, the lat- 
tice spacing A, should be several times smaller than the 
minimal De Broglie wavelength; we have found that the 
formula 



0.3 • min ■ 



1 



2 /12 



works well enough producing relative errors of order 
10- 4 . 

The fact that the discretization corresponds to the 
relative rather than absolute errors can be understood 
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as follows. The equations (|B2[) and (|B4|) may be re- 
garded as the ones describing reflection of a quantum 
particle in a kind of crystal. Indeed, the substitution 

4> k = (l- ^A^j ip k brings Eq. (|B2j) into the form 

Ti^cj) = 0, where <fi is the column composed of and 
Ti^ is a Hermitean linear operator. The probability of 
reflection in a crystal is exponentially small due to the 
same dynamical reasons as in the continuum case. One 
concludes that the finite value of A gives rise to cor- 
rections to the suppression exponent, rather than to the 



reflection probability itself, i.e. it produces relative dis- 
cretization errors. We have checked the above physical 
considerations by performing calculations on the lattices 
with different cutoffs and lattice spacings. The overall 
conclusion is that, indeed, the discretization effects re- 
sult in relative errors of order 10 -4 , while the truncation 
errors are always negligible. We kept the round-off errors 
under control by exploiting the current conservation law 
(|B3[) . which was checked to hold with precision better 
than lCT 12 . 



[1] A. M. Perelomov, V. S. Popov and M. V. Terent'ev, 

ZHETF 51, 309 (1966). [20 

V. S. Popov, V. Kuznetsov and A. M. Perelomov, ZHETF 

53, 331 (1967). [21 

[2] W. H. Miller, Adv. Chem. Phys. 25, 69 (1974). 

[3] I. Y. Kobzarev, L. B. Okun and M. B. Voloshin, Sov. J. [22 
Nucl. Phys. 20, 644 (1975) [Yad. Fiz. 20, 1229 (1974)]. 
S. R. Coleman, Phys. Rev. D 15, 2929 (1977) [Erratum- [23 
ibid. D 16, 1248 (1977)]. 

[4] M. Davis and E. Heller, J.Chem.Phys. 75, 246 (1981). [24 

[5] M. Wilkinson, Physica 21D, 341 (1986). 

S. Takada and H. Nakamura, J. Chem. Phys. 100, 98 

(1994) . [25 
S. Takada, P.N. Walker and M. Wilkinson, Phys. Rev. A 
52, 3546 (1995). 
S. Takada, J. Chem. Phys. 104, 3742 (1996). [26 

[6] T. Banks, C. M. Bender and T. T. Wu, Phys. Rev. D 8, 

3346 (1973). [27 
T. Banks and C. M. Bender, Phys. Rev. D 8, 3366 (1973). [28 
[7] V. A. Rubakov, D. T. Son and P. G. Tinyakov, Phys. 

Lett. B 287, 342 (1992). [29 
[8] A. N. Kuzne tsov and P. G. Tinyako v, Phys. Rev. D 56, 

1156 (1997) [arXiv:hep-ph/9703256j . [30 
[9] F. Bezrukov, D. Levkov, C. Rebbi, V. A. Rubakov 
and P. Tinyakov, Ph ys. Rev. D 68, 036005 (2003) 
arXiv:hep-ph/0304180 ; Phys. Lett. B 574, 75 (2003) [31 
arXiv:hep-ph/0305300 . 
[10] D. G. Levkov and S. M. Sibiryakov, Phys. Rev. D 71, [32 
025001 (2005) [arXiv:hep-th/0410198| ; JETP Lett. 81, 
53 (2005) [Pisma Zh. Eksp. Teor. Fiz. 81, 60 (2005)] [33 
[arXiv:hep-th/0412253| . 
[11] I. Affleck, Nucl. Phys. B 191, 429 (1981). 
[12] M. V. Berry and K. E. Mount, Rep. Prog. Phys. 35, 315 [34 
(1972). 

[13] S. Adachi, Ann.Phys. 195, 45 (1989). 
[14] A. Shudo, K. S. Ikeda, Phys.Rev.Lett. 76, 4151 (1996). [35 
[15] A.D. Ribeiro, M.A.M. de Aguiar, M. Baranger, Phys. 

Rev. E 69, 066204 (2004). [36 
[16] F. Parisio, M.A.M. de Aguiar, J. Phys. A: Math. Gen. 

38, 9317 (2005). 

[17] O. Bohigas, S. Tomsovic and D. Ulhno, Phys. Rep. 223, [37 
43 (1993). 

[18] E. Doron and S. D. Frischat, Phys. Rev. Lett. 75, 3661 [38 

(1995) ; 

S. D. Frischat and E. Doron, Phys. Rev. E 57, 1421 [39 
(1998). 

[19] A. Mouchet, C. Miniatura, R. Kaiser, B. Gremaud, 



D. Delande, Phys. Rev. E 64, 016221 (2001). 

A. Shudo, K.S. Ikeda, Phys. Rev. Lett. 74, 682 (1995); 

Physica D 115, 234 (1998). 

A. Shudo, Y. Ishii and K. S. Ikeda, J. Phys. A 35, L225 
(2002). 

T. Onishi, A. Shudo, K.S. Ikeda and K. Takahashi, Phys. 
Rev. E 68, 056211 (2003). 

T.Van Voorhis, E.J. Heller, Phys. Rev. A 66, 050501(R) 
(2002). 

B. Eckhardt, C. Jung J. Phys. A 19, L829 (1986); 

B. Eckhardt, Physica D 33, 89 (1988); 

P. Gaspard, S.A. Rice, J. Chem. Phys. 90, 2225 (1988). 
V. I. Zakharov, Nucl. Phys. B 353, 683 (1991); Phys. 
Rev. Lett. 67, 3650 (1991). 

G. Veneziano, Mod. Phys. Lett. A 7, 1661 (1992). 

M. Maggiore and M. A. Shifman, Phys. Rev. D 46, 3550 

(1992). 

M. B. Voloshin, Phys. Rev. D 49, 2014 (1994). 

V. A. Rubakov and D. T. So n, Nucl. Phys. B 424, 55 

(1994) [arXiv:hep-ph/9401257| . 

C. S. Drew, S. C. Creagh and R. H. Tew, Phys. Rev. A 
72, 062501 (2005). 

F. Bezrukov and D. Levkov, arXiv:quant-ph/0301022 J. 
Exp. Theor. Phy s. 98, 820 (2004) [Zh. Ek sp. Teor. Fiz. 
125, 938 (2004)] [arXiv:quant-ph/0312144] . 
K. Takahashi and K.S. Ikeda, J. Phys. A 36, 7953 (2003); 
Europhys. Lett. 71, 193 (2005). 

F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 
2212 (1984). 

G. F. Bonini, A. G. Cohen, C. Rebbi and V. A. Ru bakov, 
Phys. Rev. D 60, 076004 (1999) [arXiv:hep-ph/9901226l ; 
|arXiv:quant^ph/ 9901062 

D. G. Levkov, A. G. Panin and S. M. Sibiryakov, 
"Tunneling via unstable semiclassical solutions," 

larXiv:0707.0433l [quant- ph] . 

D. G. Levkov, A. G. Panin and S. M. Sibiryakov, Phys. 
Rev. A 76, 032114 (2007) [arXiv:0704.0409] . 
W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flan- 
nery, "Numerical recipes in C: the art of scientific com- 
puting", Cambridge University Press, 1992. 
A.L. Xavier Jr., M.A.M. de Aguiar, Ann. Phys 252, 458 
(1996). 

D. Levkov, C. Rebbi and V. A. Rubakov, Phys. Rev. D 
66, 083516 (2002) [arXiv:gr-qc/02060"28] . 
\http://solver. inr. ac.ru 



